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SUMMARY 


A procedure is presented for effective consideration of viscous effects in 
computational development of high Reynolds number flows. The procedure is based 
on the interpretation of the Navier-Stokes equations as vorticity transport 
equations. The physics of the flow has been represented in a form suitable for 
numerical analysis. Lighthill's concept for flow development for computational 
purposes has been adapted. The vorticity transport equations were cast in a 
form convenient for computation. A statement for these equations was written 
using the method of weighted residuals and applying the Galerkin criterion. An 
integral representation of the induced velocity was applied on the basis of the 
Biot-Savart law. The numerical procedure developed by Hess and Smith and then 
refined by Argyris and Scharpf was used in computing vortex sheets over curved 
wing surfaces. Distribution of new vorticity, produced at wing surfaces over 
small computational time intervals, was assumed to be confined to a thin region 
around the wing surfaces. This report presents thus the foundation and compo- 
nents of the numerical method. A following report will present the details of 
the implementation of the computational procedure. 



INTRODUCTION 


Theoretical treatment of unsteady viscous flow effects and extension of the 
Reynolds number range are of primary concern as well as having great urgency 
from an application point of view (ref. 1). For higher Reynolds number flows, 
the effect of viscosity becomes increasingly confined to narrow regions of the 
wing surface and laminar flows break down and become turbulent (ref. 2) . The 
treatment of viscous effects and high Reynolds number flows is still not pre- 
sented well in theoretical analyses (ref. 1). 

The Navier-Stokes equations, as the mathematical model of unsteady viscous 
flow, contain all features needed for the solution of the problem. Advances, 
however, in establishing practical, meaningful relations between vorticity de- 
velopment and diffusion have been slow. Ingenious schemes have been applied 
in seeking to obtain solutions of the Navier-Stokes equations. Solutions of 
this kind have thus far been confined to flows of relatively low Reynolds num- 
bers, and they are still inadequate for actual aeronautical applications (ref. 

1) . By applying suitable discretization and computational procedures, attempts 
have been made to extend the solutions to higher Reynolds numbers (ref. 3) . The 
obtained results with viscous flows of moderately high Reynolds numbers have 
been achieved through a long struggle with mathematics, fluid dynamics, and te- 
dious numerical analyses. For instance, presented data and computer plots, al- 
though giving a blurred picture of viscous flow patterns, have been of genuine 
interest (ref. 4) . 

Since the finite-element methods are well suited for theoretical treatment 
of unsteady viscous flows past bodies of complex geometry, and since config- 
urations of finite-element meshes can be suitably clustered in regions of sig- 
nificant changes in the flow patterns, there has been a vivid interest and much 
work done in this area (ref. 5) . There is expectation that finite-element meth- 
ods will provide a solid numerical foundation for treatment of high Reynolds 
number flows. The purpose of this report has been the presentation of such a 
foundation. 

It appears that four different formulations have been commonly applied in 
conjunction with discretization methods: the formulations with respect to stream 
function, stream function and vorticity, velocity and pressure (ref. 6), and the 
Lighthill vorticity-velocity concept of flow development for computational pur- 
poses (ref. 7). Interpretations of the Lighthill concept, using the finite-dif- 
ference method, have been presented by Payne (ref. 8) and by Schmall and Kinney 
(ref. 9), as well as in conjunction with an integro-differential formulation by 
Wu (refs. 10 and 11). 

The authors are deeply grateful to Dr. G.M.L. Gladwell, Professor, Univer- 
sity of Waterloo, Waterloo, Ontario, Canada, and to Dr. Ch. Hirsch, Professor, 
University of Brussels, Brussels, Belgium, for the kindly provided assistance 
and specialized expertise for the research toward this report. Acknowledged 
gratefully is also the support from the NATO grant No: 1264 through which it 

was possible for Professor Hirsch to participate in the research. 
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SYMBOLS 


A 

dA 

c 


g 

L. 

J 


£ 1 ,£ 2 ,£ 3 

M k 


m 


N. 


1 


n 

fluid region, m 

2 

fluid element containing vorticity w(t+6t), m 
wing boundary 

number of integration points 
Lagrangian interpolation polynomials 

area coordinates 

weighting functions 
number of element nodes 
nodal shape functions 


n unit vector normal to the wing surface 

parameters 


dS 

t 

t 

fit 

U 


v(t) 

v(t+6t) 


v ft+6t) 
-Y 

v 
— 0 ) 


V 

—to 


* 


position vector relative to dA or dS respectively; also a distance, 
element of a vortex sheet 
time, sec 

unit vector tangent to the wing surface 
subsequent time-step, sec 

uniform onset velocity of the fluid in region A; constant in time, 
m/sec 

nodal values of velocity, m/sec 
no-slip velocity, m/sec 

velocity which does not satisfy the no-slip condition at the begin- 
ning of iteration, m/sec 

induced velocity due to vortex sheet, m/sec 
induced velocity at a point in the fluid, m/sec 
additional induced velocity, m/sec 


W a function 

w^ weights 


x,y,z Cartesian coordinates 

nodal values of vortex sheet intensity, m/sec 

yCt+fit) vortex sheet, m/sec 



V 


2 

kinematic viscosity, m /sec 
£,n element local coordinates 

3 

p uniform density, kg/m 

ok vorticity at a finite-element node, 1/sec 

oj (t) vorticity, 1/sec 

o)* additional vorticity during time <5t, 1/sec 

oj* undetermined parameters 

Matrix Symbols 

[A] convection coefficient matrix 

{B} finite-element vector containing wing boundary conditions 

[C] vortex sheet coefficient matrix 

{F} a vector 

{y} unknown vector 

[J] Jacobian 

[m] mass matrix 

[S] diffusion matrix 

Subscripts 

e over a finite-element 

f.e. finite-element 

i element node 

x,y partial derivatives with respect to x and y respectively 
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THE MATHEMATICAL ANALYSIS 


In unsteady viscous flows around wing boundaries, there exists a very nar- 
row region in which gradients of flow field variables and influence of viscous 
forces are very large. At higher Reynolds number flows, and flows around os- 
cillating wings in particular, the situation is even more complicated. The 
proper solution requires consideration of the entire flow field around the wing 
There is a need of a reasonably efficient and accurate method for treatment of 
such flows. In principle, the time -dependent Navier-Stokes equations, as the 
mathematical model, contain all the features needed for the solution of the 
laminar and turbulent aspects of the flow. The interpretation of these equa- 
tions as vorticity transport equations extends Helmholtz's equations to vis- 
cous flows (refs. 12 and 13). 

Numerical methods of solution appear to provide the only hope of producing 
quantitative results. The procedure described below offers the foundation of 
such an approach. For reliable results about the viscous region around the 
wing, the Lighthill concept of flow development for computational purposes 
(ref. 7) has been adapted here. It is thus important to provide for the prop- 
er inclusion of vorticity in the numerical analysis of the field variables of 
the complete flow. 


Summary of the Computational Procedure 

The vorticity transport equations were cast in a form convenient for compu 
tation (refs. 14 and 15). The computation is considered unsteady in the sense 
of being suitable for time-stepped solutions. A statement for these equations 
was written using the method of weighted residuals and applying the Galerkin 
criterion (see appendix A). From the discretized form of the equations, finite 
element equations were assembled and summed over all elements forming thus the 
global system of equations for the complete flow region. The discretization in 
finite-element form of the vorticity transport equations resulted in a system 
of differential equations with respect to time. These equations can be solved 
with respect to the time variable using a finite difference technique (ref. 5) . 
Finite-elements in the fluid region contain vorticity, since vorticity will be 
found only in that portion of the fluid region which is discretized by finite- 
elements . 

An integral representation of the induced velocity at finite-element nodes 
due to the vorticity distribution was applied on the basis of the Biot-Savart 
law. Details are shown in appendix B. The numerical procedure developed by 
Hess and Smith (ref. 16) and then refined by Argyris and Scharpf (ref. 17) was 
used in computing vortex sheets over the curved wing surface (see appendix C) . 
The refinement introduced by Argyris and Scharpf involved curvilinear elements, 
in place of the plane elements used by Hess and Smith (see appendix D) . A 
Hermitian interpolation model was used to approximate the geometry of curvilin- 
ear elements between nodes defining the elements. A Lagrangian interpolation 
model was used to describe vortex sheet intensity over finite -elements . Two 
types of finite-elements were applied in representing the flow region under in- 
vestigation: isoparametric triangular and isoparametric quadrilateral finite- 

elements. Surrounding the curved wing surface are quadrilateral finite- 
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elements, with the remainder of the fluid region divided into triangular fi- 
nite-elements of varying sizes, becoming smaller and smaller near the wing sur- 
face. Quadratic shape-functions were utilized for both isoparametric triangular 
and isoparametric quadrilateral finite-elements. Details are shown in appendix 
E. 

Vorticity transport equations - In three-dimensional form the vorticity trans- 
port equations can be written in fixed coordinates (ref. 13) as 

— = (a)/V)v - (v»V)(£ + vV 2 io (1) 

du_ 

In the case of two-dimensional flow Qj*grad)v = 0 and -— = vV 2 co, i.e., the 

total variation of vorticity is equal to the rate of dissipation of vorticity 
through friction (ref. 12) . The two-dimensional form of equation (1) is 

— = - (v-V)io + vV 2 (o (2) 


whereby only convection and viscous conduction take place (ref. 15) . Equation 
(2) is discretized in space using the finite-element method. Over each finite- 
element, the vorticity and velocity components are approximated by their nodal 


values. 


to . 

J 


u. 


and v., and by the nodal shape functions 


N. 

J 


Then 


m 

<*>(£, n) = E 
j=l 

Nj (5,n)a)j 

= (N} t {a J } e 

(3) 

m 

u(?,n) = E 
i=i 

N. ce,n)u 

= {N} t {u) e 

(4) 

m 

v(S,n) = E 
j=l 

N. (S,n)v. 

= (Nr{v} e 

(5) 


The finite-element is defined here by m nodes. The element shape functions 
are expressed in terms of the local coordinates (£,n), which range between 
unity and zero within the finite-element. The coordinate r) will be chosen 
to be normal to the wing boundary. The weighted residual statement for the 
vorticity transport equations is written as 

¥s*» ^ * 37 (v “ ) - v ( * w J 1 ’ 0 (6) 

In equation (6) the weighting functions of the method of weighted residuals have 
been set equal to the finite-element shape functions, N^, in accordance with 

the Galerkin criterion. Applying an integration by parts to the second-order 
viscous terms, the expression becomes 
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( 7 ) 



tvfr + h w 


* (VU)} + V 


BN. a BN. a 

r 1 jko + 1_ _3<d i 

i 3x 3x 3y 3y 


dA 



0 


The expression in equation (7) is discretized in finite-element form. Over each 
finite-element, the vorticity transport equations can be written in matrix form 
as 


[m] 6 {~} e = [A] e {u)} e + [S] e {m} e + {B} e 


( 8 ) 


where the element matrices and the vector, after substitution of equations (3), 
(4), and (5) into equation (7) are defined as follows: 


[m] e 



{N} {N}*' 


dA 


( 9 ) 


[A] 


e 



[{N}{u} t {N}{N l 1 + (N}{v} t {N}{N l 1 ] dA 
x y 


( 10 ) 


[S] & = V J A [{N x HN x }t + {N y HN y }t J dA (11) 

Here equations (9-11) are respectively the element mass matrix and the convec- 
tion and diffusion matrices. Matrix [A] is an unsymmetric matrix, representing 
the nonlinear convective terms in the vorticity transport equations. It is a 
function of the no-slip velocity field v_(t) corresponding to the no-slip con- 
dition. Partial derivatives with respect to x and y are indicated by sub- 

scripts x or y respectively. The finite-element vector {B} , containing 
the wing boundary conditions, is written as 


(B) 


e 



(N) (|^}dc = [B] e {m} 6 


( 12 ) 


It results from the integration by parts. Here 

[B] 6 = v (NHlO 1 dc 

The element matrix [B] is constant in time and is computed only once for each 
element and then stored. The contour integral is not evaluated on the outer 
boundary since the region covered by finite-elements contains all the vorticity 
in the flow and therefore m is identically zero at the outer boundary. 

The numerical integration used to evaluate the integrals of equations 
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(9-12) is carried out in the local coordinate system of the finite-element, 
where 


dx dy = det. [J] d£ dn (13) 

is used to transform the region with respect to which the integration is carried 
out. The vorticity transport equation, in the form of equation (8), is then 
solved in time over the subsequent time-step 6t using a finite-difference 
technique, to obtain a new distribution of the vorticity at time t+6t, which, 
in general, is different from that vorticity distribution which is known at time 
t. 


THE DEVELOPMENT OF THE FLOW 


The manner of progressing from the velocity and vorticity distribution at 
an initial time instant t to the vorticity and velocity distribution of a sub- 
sequent time instant (t+St) will now be outlined. Consider, however, that 
some time t has elapsed since the fluid was set in motion past the wing. 

The progressing in time - At time instant t, the flow variables, v;(t) and 
co(t), are assumed to be known and that the no-slip velocity condition is sat- 
isfied on the wing boundary c. That is, in solving equation (8) the no-slip 
velocity field v_(t) has to be used. Since vorticity is a measure of the an- 
gular momentum of the fluid particles in motion, in the two-dimensional case 
the vorticity to has a single component, perpendicular to the plane of motion 
x and y. 

During the subsequent time interval 6t, the flow variables v_(t) and 
co(t) are advanced to their new values at t+<5t. Through convection and diffu- 
sion, the known vorticity <o(t) is redistributed, while during the same time 
interval, St, new vorticity is diffusing into the fluid from the wing bound- 
ary c in such a manner that the no-slip condition is satisfied. 


The Biot-Savart law - Consider first the redistribution of the known vorticity 
m(t) through convection and diffusion. By solving equation (8), the redistri- 
bution of to(t) at t+St is given as co(t+St). Then the induced velocity, 
which coexists along with the vorticity io(t+6t), is given by the Biot-Savart 
law as 


v 

— u 


(t+6t) 


1 _ 

277 



w(t+6t) x r 



dA 


(14) 


where the integration is over the whole region of the vorticity distribution. 
Here t is a position vector relative to the fluid element dA containing vor- 
ticity aj(t+6t). 

Conversely, if v^(t+St) is the induced velocity at a point in the fluid, 
the vorticity at that point is by definition 

w(t+6t) = x v (t+6t) (15) 
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Here, v^ (t+fit) is computed by numerically integrating equation (14) over each 
finite-element containing vorticity and summing their individual contributions. 

Development of the vorticity - In general, the induced velocity v^ (t+fit), 

resulting from the vorticity as obtained from equation (8) , has a non-zero slip 
component at the wing boundary. Therefore, the velocity v(t+6t) does not, in 
general, satisfy the no-slip boundary condition on the winjf surface at the be- 
ginning of the iteration process. The difference from the no-slip condition 
can be attributed to some additional vorticity <d* ( t+fit) , produced at the wing 
surface over the time step fit but which was not included when the induced ve- 
locity was computed, i.e., 

io(t+fit) + to*(t+5t) = V_ x v^ (t+fit) + V_ x v^*(t+5t) 

Accordingly, any existing velocity field can be regarded as containing two in- 
duced velocity vectors, v^(t+St) and v^*(t+fit). The complete velocity field 

at t+fit, at any point in the fluid, can be defined (refs. 7, 18, and 19) as 

v(t+6t) = v^t+fit) + (t+fit) + (16) 

The velocity ;v (t+fit) , defined by equation (16), should satisfy the no-slip 
condition at the wing surface. The induced velocity v^ (t+fit), obtained from 

viscous considerations, is unequal to zero at the wing surface. It is impor- 
tant to consider that in the framework of the concept of computational flow de- 
velopment, used in obtaining the solution for the flow by numerical means, we 
simulate the flow of vorticity in the given unsteady viscous flow by introduc- 
ing elements of irrotational. flow. The additional induced velocity v^* (t+fit), 

although calculated from irrotational flow considerations, is a viscous effect, 
and it is rotational by definition. It is also unequal to zero at the wing sur- 
face. 


Consider now the resulting tangential velocity at the wing surface (ref. 7), 
before the new vorticity w* (t+fit) was accounted for. Since it may not satis- 
fy the no-slip condition there, it appears that just enough new vorticity 
w* (t+fit) must have been diffused at the wing surface during the time interval 
fit, that the sum of the tangential components of the induced velocities due to 
u (t+fit) and w* (t+fit) and the onset velocity will lead to zero slip. Corre- 
spondingly, equation (16) may be rewritten for the tangential components of the 
involved velocities on the wing surface. Upon rearranging the resulting equa- 
tion one obtains 


t • v . (t+fit) = - t . [v (t+fit) + U 1 (17) 

In sum, given the vorticity to (t+fit) from the solution of equation (8), 
the induced velocity v^ (t+fit) and the velocity ;v (t+fit) can be determined 

from equation (16). The velocity v_(t+fit), however, may not satisfy the no- 
slip condition on the wing surface c. Therefore, an additional distribution 
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qf vorticity w*(t+6t) must be combined with the existing vorticity distribu- 
tion w(t+6t). The induced velocity v (t+6t) can be regarded also as a per- 
turbation velocity (refs. 18 and 19). 


Intensity of the new vorticity - As iteration for time-step t+6t begins, the 
tangential components of velocity v(t+6t) at the wing surface may be unequal 
to zero. There exists then a discontinuity in the tangential velocity compon- 
ent at the wing surface, equal to the difference between the no-slip velocity 
and the velocity component t_*v_(t+5t).. By definition, a surface across which 
tangential velocity changes abruptly is a vortex sheet (ref. 18). The concept 
of a vortex sheet is introduced now as an approximation, allowing the calcula- 
tion of the additional vorticity oj*. Following reference 18, with the in- 
duced velocity v_ (t+6t) due to the vortex sheet y(t+<5t), the velocity at a 

point P in the fluid can be written now as 


v^(t+6t) = v^(t+<5t) + v** (t+St) + U 


(18) 


Following the derivation of the vortex sheet (ref. 18) , the induced veloc- 
ity due to the vortex sheet at a point I, just outside the wing surface, and 
the velocity at a point 0, just inside the wing surface, may be written as 


v I (t+6t) = v*(t+6t) + y yft+iSl;) x n (19) 

y y ^ 

and 

v°(t+6t) = v*(t+6t) - y- Y_(t+<St) x n (20) 

respectively. Since a straight element of vortex sheet dS cannot induce ve- 
locity at its centroid (ref. 18), the velocity v^(t+6t) is the velocity in- 
duced at its centroid by the rest of the vortex sheet, excluding the infinites- 
imal element dS. Moreover, the second terms on the right-hand side of equa- 
tions (19) and (20) are the velocities induced at points I and 0 respective- 
ly by the element dS (ref. 18) . 

If a general point P is now identified as point 0, equation (20) may 
be substituted into equation (18) as 

v^(t+6t) = v^(t+6t) + v*(t+6t) - y y (t+iSt) x n + (21) 

(i) y “ 

Here ii is the unit vector normal to the wing surface. Furthermore, by the 
above definition of the vortex sheet, the velocity at point 0 is set equal to 
the no-slip velocity, 

v^(t+6t) = 0 = v^(t+6t) + vj£(t+6t) - j y(t+6t) x n + (22) 

Taking the tangent component of equation (22) and rearranging 


10 



(23) 


t_ • v^(t+5t) - y Y(t +| St) = - t_ • (U^ + v^(t+<5t)) 


The velocity v*(t+6t) has been given (ref. 18) as 


v*(t+6t) 


= f s 


Y (t+6t) x r 


dS 


(24) 


which is the two-dimensional form of the Biot-Savart induction law for vortex 
sheets. Substituting equation (24) into equation (23) 


t • 



Y (t+6t) x r 1 „ 

v-ri dS - ^ Y (t+St) = - t * [^(t+et) + Uj 


(25) 


Equation (25) is solved to obtain the vortex sheet. The integration is carried 
out over the wing boundary (ref. 18). The amount of new vorticity aj*(t+6t), 
diffusing from the wing surface during the time interval 6t, is found using 
the vortex sheet y(t + dt). A sufficient amount of new vorticity is created 
close to the wing boundary such that when combined with the existing vorticity 
oi(t+6t), the no-slip condition is insured. This condition, explained by Light- 
hill (ref. 7), insuring that the velocity normal to the airfoil surface will be 
zero, is expressed analytically here as follows 


t+St 



Y (t+St) 


(26) 


Nearly similar interpretations of this condition have been expressed analyti- 
cally in references (9-11) . Since diffusion of vorticity from a vortex sheet 
can be considered as a time-dependent problem, the distribution of new vortic- 
ity w*(t+6t) in the fluid could have been carried out using the finite-ele- 
ment vorticity transport program described previously. Since the process is 
iterative, it is carried out until the velocities around the wing are equal to 
zero for a given numerical accuracy. However, in this report an alternative 
approach will be used. 


When it can be assumed that diffusion dominates only the vorticity 
w*(t+<St) over convection near the wing surface and not the vorticity a)(t+6t) 
and provided that the time-step is sufficiently small (ref. 7), the new vortic- 
ity co*(t+6t) becomes (ref. 15) 

-(n-r) 2 , 


u)*(t+6t) = 


Y (t+6t) 


4vSt 


{irv<5t} 


(27) 


Thus, included in equation (27) is the idea that at time t+6t, the vorticity 
co*(t+6t) on the left-side of the equation and the vortex sheet Y(t + <$t), on the 
right-hand side, occur simultaneously. The range of validity of equation (27) 
is confined to the region near the surface where viscous effects are much 
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larger than the vorticity convection (i.e., a computational thickness, y, 
much smaller than the boundary layer thickness, <S) and for times when the time 
derivative of diffusion is much larger than the convection of vorticity (i.e., 
when U^St is much smaller than the chord multiplied by 6/y) . 

DETAILS OF THE SOLUTION PROCEDURE 
The Vorticity Transport Equations 

The details of the implementation of the solution of the vorticity trans- 
port equations, in the form of the discretized equation (8), can be summarized 
as follows: 

a) the global mass matrix [m] is constant in time for a given finite- 
element mesh. It is assembled and stored in a symmetric band-storage 
mode at the beginning of the numerical procedure. This matrix is 
computed using numerical integration and then decomposed into equiva- 
lent lower and upper triangular matrices using Cholesky decomposition. 

b) the time rate of change of vorticity vector {—} is obtained using 

o X, 

a finite-difference scheme; the initial procedure involved the Euler 
method with later use of the trapezoidal method and the fractional- 
step method. 

c) the finite-element convection coefficient matrix [A] is dependent on 
time due to the velocity terms and must be evaluated using numerical 
integration each time the vorticity transport equations are solved. 

d) the finite-element diffusion coefficient matrix [S] is constant in 
time. It is computed using numerical integration at the beginning 
of the computational procedure and then stored. 

e) the matrix [B] is assembled for the wing boundary finite-elements 
only and is evaluated using numerical integration; it is also con- 
stant in time and stored at the beginning of the computational pro- 
cedure . 


Computation of the Additional Vorticity 

The additional vorticity can be distributed using the finite-element vor- 
ticity transport program or the solution due to Lighthill’s approximation, 
equation (27) . The first approach can easily be implemented in the computer 
code. Hence, in what follows in this report, only the way of applying equa- 
tion (27) will be described. Assumptions can be made that the distribution of 
new vorticity w*(t+6t), produced at the wing surface over a small time-inter- 
val 6t, is confined to a thin region around the wing, and that it is dominat- 
ed by diffusion. 

The distribution of new vorticity determined using Lighthill's approxima- 
tion, as expressed by equation (27), is proportional to the vortex sheet. A 
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vortex sheet may be computed from equation (25) . The second term on the left- 
hand side of equation (25) indicates that the singularity has been properly con- 
sidered. The form of the right-hand sides of both equations (17) and (25) is of 
interest. Accordingly, the right-hand side of equation (25) is the negative of 
the resulting slip velocity v. (t+6t) at a point of the velocity field v(t+5t). 

This slip velocity has been defined (ref. 18) as a discontinuity in the tangen- 
tial component of the velocity at the wing surface. 

Equation (25) is solved numerically at a finite number of points on the 
wing surface. The vortex sheet y(t+St) is described by these finite number 
of points and suitable interpolation functions (refs. 16 and 17). The inte- 
gral in equation (25) is evaluated using a Gauss -Legendre quadrature and as- 
sembled for all nodes on the wing surface. The resulting system of equations 
is solved for vortex sheet y(t+6t) using a LU-decomposition. Once the vor- 
tex sheet has been computed, the distribution of additional vorticity is ob- 
tained. As it will be shown in the follow-up of this report, the actual com- 
putation of the additional vorticity may be confined to only the nodes of the 
first few layers of finite-elements around the wing. 

At time t+<5t, the new vorticity u)*(t+<5t) is combined with the pre- 
viously computed vorticity u>(t+6t), with the result becoming a new estimate 
of the existing vorticity w(t+<5t) at time t+St or 

w(t+5t) u>(t+6t) + io*(t+St) (28) 

Equation (28) is evaluated at all finite-element nodes where the new vorticity 
oo*(t+6t) was computed. 

The new estimate for the existing vorticity ui(t+6t) is then used to re- 
calculate the velocity v(t+6t) from equation (16) . This velocity should sat- 
isfy the no-slip condition at the wing boundary. However, it may not, since 
the additional vorticity distribution ui*(t+6t) is only an estimate based on 
the calculated value of the vortex sheet y(t+<5t). Thus, if it does not, 
equation (25) can be solved again for a new vortex sheet y(t+St), with the 
subsequent solutions of equations (27) and (28) giving still another estimate 
of the existing vorticity io(t+<5t) at time t+6t. Equations (25), (27), and 
(28) may be iterated in this manner until a sufficient amount of new vorticity 
w*(t+St) is combined with the existing vorticity co(t+St) to enforce the no- 
slip condition at the wing boundary. Once the no-slip condition on the wing 
boundary is satisfied, the solution is complete and the computational proced- 
ure may be extended over the next small time-step, t+6t; starting with the 
solution of the vorticity transport equations. Figure (1) illustrates the 
main sequences of computations in the computational procedure. 

CONCLUDING REMARKS AND RECOMMENDATIONS 

The presented procedure offers an approach toward effective representa- 
tion of viscous effects in the computational development of high Reynolds num- 
ber flows. Efforts were made to represent the physics of the flow in a form 
suitable for numerical analysis, to contain the features needed for the 
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solution, and to optimize codes to reduce computer costs. Since finite-ele- 
ment methods are considered as computationally effective, it is expected that 
they will facilitate extending not only the range of Reynolds numbers of prac- 
tical applications, but they hold promise for good results from three-dimen- 
sional analyses. 

The numerical difficulties increase as the Reynolds numbers increase; an- 
alytical difficulties adversely affecting accuracy. Computer runs were made 
to test the range of application of the program and the accuracy of the results 
for higher Reynolds numbers. These tests have required a considerable amount 
of computer execution times. Portions of the program need significant optim- 
ization in order to reduce computer cost. The computation, for instance, of 
the induced velocity from the Biot-Savart law still involves a large portion of 
the computer time for a time-step. Significant instabilities occur in the nu- 
merical integration which are related to the size of the time-step of integra- 
tion. Integration and computational techniques have been employed in the pro- 
gram with the objective of increasing the integration time-step. It is expected 
that the implementation of the fractional -step method will further reduce com- 
puter execution times. The numerical treatment of the wake region requires the 
application of the correct boundary conditions at the downstream end of the fi- 
nite-element gridwork. It is thus necessary that these conditions account for 
the propagation of vorticity over the finite-element gridwork. 

Particularly promising, for instance, is the Boundary-Element Concept which 
was developed just recently as a refinement in the field of discretization tech- 
niques. It has also been necessary to determine how the numerical results sim- 
ulate the actual flow at comparable Reynolds numbers. There is available only 
a very limited amount of suitable experimental results for comparisons. How- 
ever, the prospects for such comparisons have improved since interest in suit- 
able experimental research has also increased. 
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Figure 1. The main sequence of computations in the computational procedure. 

15 

















Figure 1, Concluded. 








APPENDIX A 


The Vorticity Transport Equations 

The vorticity transport equations are discretized in space employing the 
finite-element method. Element equations are defined from a functional and 
then assembled to form a matrix-system of first-order differential equations. 
This system governs the fluid region. Using the method of weighted residuals, 
a statement for the vorticity transport equations can be written as 


L 


... ( 3(i) 3&) , 3u) 

W { 3t * U Df * V 3? 


f 3 2 u) 

V C 3^ + 


3fw 

3y 


T ) ) dA = 0 


(Al) 


The function W is approximated as 


m 

w = z 

k=l 




(A2) 


k k 

where the are arbitrary parameters and the M are weighting functions. 

On substitution of equation (A2) into equation (Al) , the statement becomes 



m 

£ 

k=l 


„k k 

M u>. 


1 3d) 
* 1 3t 




Since the parameters 


0). 


are arbitrary 



+ 


3 2 <d 

ay 2 " 


} dA 


0 (A3) 




k=l 



3co 

¥x 


+ 





} dA 


0 (A4) 


In equation (A4) , the highest order derivatives are the second-order terms 
multiplied by v, the viscous terms. Taking only the viscous terms from equa- 
tion (A4) and applying a transformation by integration by parts (ref. 5) re- 
sults in 
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On substitution of equation (A5) back into equation (A4) , the statement for 
the vorticity transport equations becomes 


/ 


E [ M k { 9u) 

K=1 


3 co 

at + u ax 


3w 

V F } + 


v { 


3M k 

ax ax 


3M k 3w 

ay ay 


} ] dA 
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k=l 


M 1 dc 
3n 


(A6) 
= 0 


where the highest order derivatives are now of the first-order. Applying the 

v 

Galerkin process to equation (A6) , the weighting functions M are set equal 
to the finite-element shape functions N^. 

The first space derivatives of the vorticity are determined from equations 
(E13) as follows — 


ao>(g,n) 

3x 


m aisLC^n) 


E 

i=l 


8x 


w . 

1 


(N } t {w} 6 
X 


(A7) 


and 


a^(g,n) = 

ay 


m 3N. (5 ,t0 

i?! -5-“i 


{N } t {w} 6 


CA8) 


Similarly, the vorticity time-derivative is written as 


- V N f£ nl — = {N}* 

at . . i u,nJ at 11 at 

1=1 


(A9) 


Substitution of equations (A7-A9) and equations (E13-E15) into equation (A6) 
results in a matrix-system of differential equations over one finite-element as 


[m] 


e 8{w} 


3t 


[A] e {w} 6 + [S] 6 {w} e + {B} 6 


CA10) 


The element mass matrix [m] is given as 


[m] 


L 


{NHN} dA 


(All) 


and the element convection and diffusion matrices are given as 
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[A]' 


7 [{N}{u} t {N}{N }* + {N}{v} t {N}{N ] dA (A12) 

A X y 


[S] c = - v 


/ [{N X HN x }t + {N y HN y }t 3 dA 


(A13) 


respectively. 

conditions 


The element vector {B} contains the wing vorticity boundary 


{B} = v 


/ 


{N} dc 

dn 


(A14) 


The integration of equations (A11-A14) being very difficult, it is therefore 
replaced by a numerical integration procedure as follows 



f(x,y) 


dA 


g 

r £ <w wC W 

k=l 


(A15) 


Here the terms W (S] C > , \) are t ^ ie integration weights corresponding to the total 

number of g integration points £ k ,n k . The integration over a finite-element 

is carried out in the local coordinate system of the finite-element. To trans- 
form a region to the one in which the integration is carried out, i.e., from the 
x,y plane to the £,n plane 

dA = dxdy = det[J] d£dr| (A16) 

The finite-element matrices of equations (A11-A14) are thus numerically evalu- 
ated as 

g 

[m] e = £ {NC5 k ,n k )HN(5 k ,n k )} t det. [J(5 k »\)] w(5 k ,\) (A17) 

k=l 


[A] 6 =-E {{N(5 k ,n k )Hu} t {Na k ,n k )} t {N x C5 k ,n k )} t 

k=l 

(A18) 

+ {N(C k ,n k )Hv} t {N(c k ,n k )} t {N y C5 k> n k )} t } det. w (5 k >\) 
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e g 

[ 5 ] = -vE 

k=l 


CA19) 


« N x C? k ,n k )}{N x ( 5k ,n k )> t 


+ ^y^k’V^y^k'V* } [J(? k’V ] wCC k , V 


and 

g 

[B] 6 = v E {N(C k ,-l)}{N n CC k ,-l)} t [J(5 k ,-1)] w(C k ) {w} CA20) 

k — X 

respectively. In equation (A20) the integration is over the side of the finite- 
element in contact with the wing surface. For the quadrilateral finite-element 
this side corresponds to n = -1 (see figure E.2). The derivatives in the 
Cartesian coordinate system of the element shape functions are obtained in 
terms of their local derivative with equation (E5) . 

The element equations are now assembled to form a global matrix-system of 
equations, which governs the fluid region as 

[m] = [A] (a)} + [S] (a)} + [B] {oj} (A21) 

Equation (A21) is then solved in time using a finite-difference technique to ob- 
tain a new distribution of vorticity. 
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APPENDIX B 


The Velocity Field 

The velocity field v(t+6t) of the fluid region is defined in general as 

v(t+6t) = v^t+St) + (Bl) 

The induced velocity v (t+<5t) is given by the Biot-Savart law. It is due to 

the vorticity distribution us (t+St) in the fluid region. The velocity U 
is the onset velocity. It is constant in space and time. Defined by their val- 
ues at all finite-element nodes in the fluid region, the velocity components at 
a finite-element node i are given as 


u*(t+5t) = u^ + u^(t+5t) 


v (t+St) = v + v (t+5t) 

V co us v 


(B2) 


Each finite-element will contribute toward the induced velocity at a finite- 
element node i due to the vorticity distribution over it. For finite-element 
A., the contributions to the induced velocity components are obtained as 


and 


with 


u^t+St) 

US 



Wj (t+St) (y^y^) 
]r~.[ 2 

'-lj 1 




v^ (t+St) 

US 



us. (t+St) (x.-x.) 

3 1 r 




r. . 
-13 


(x.-x.) 2 + {y.-y -) 2 

i y J x y 


(B3) 


(B4) 


(B5) 


The individual velocity components are summed over all finite-elements (f. e.) 
to give the complete induced velocity at node i as 


v*(t+St) 


f.e. . f.e. 

22 u^(t+st)i^ + 22 v^(t+st)j_ 


(B6) 
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The integrals of equations (B3) and (B4) are integrated numerically, using 
equations (A15) and (A16), as follows 


/ir k=l 


0) 

and 


i i. - y(5 k >\)> 


l^ik' 


, t g w(€.,n ) {x - x(C v ,n v )> 

v£(t+«t) = £ L_JL i~ k k 


2tt 


k=l 


la , 12 


det. t J C5 k ,n k )] w (5 k >\) CB7) 


det. [ J C5 k ,n k 3] w (S k >\) (B8) 


-ik 1 


The vorticity and the coordinates x (£ k >\) and yC? k > r l k ) are ob- 

tained from equations (El), (E2) , and (E13) . For the quadrilateral finite-ele- 
ment and the triangular finite-element, the terms of the Jacobian [J] are 
computed from equations (E7) and (Ell) respectively. 
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APPENDIX C 


The Vortex Sheet 

The vortex sheet is determined using a numerical technique developed by Hess 
and Smith (ref. 16) and which was refined by Argyris and Scharpf (ref. 17) . 
Collecting the nodal values y of the vortex sheet intensity into the vector 
{y}, equation (25) may be written in matrix form as 

[C] {Y> = {F} (Cl) 

Solving equation (Cl) for the unknown vector {y} results in 

{y} = [C]" 1 {F> (C2) 


The manner in which the terms of the coefficient matrix [C] are determined is 
now summarized. Equation (25) is applied to one of the curvilinear element- 
nodes and is then integrated around the wing-contour c. The integration around 
the contour c is accomplished by integrating over segments of contour c, de- 
fined by the curvilinear elements and then summing their contributions. For 
node i and curvilinear element c^ the terms of matrix [C] are given as 




L. (x. ,y. ) 


t. 
— i 



where the components of the vector J3 are given as 


CC3) 


and 


S 

x 


S 

y 



L j( x k ,y k ) c yi - 


y k> 




L j ( v y k > < x i - V 




(C4) 


(C5) 


Here the parameters 


L. 

J 


are the Lagrangian interpolation polynomials of equa- 


tions (D8) . They define a vortex sheet of unit intensity over curvilinear ele- 
ment c^. The position vector r^ is given as 


r . 


— ik 1 



x k ) 2 ♦ C/i-y ,) 2 


(C6) 
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The terms of the vector {F}, on the right-hand side of equation (Cl), are 
determined as 


F. 

1 


- t. 
—l 


(U 


+ v ) . 
—or 1 


(C7) 


The induced velocity v^ is obtained from equations (B7) and (B8) . Equations 
(C4) and (C5) are integrated numerically as follows 


and 
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x 


g 

£ 

k=l 


- * e k> } 





1 _ 

2 w, 
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+ 


g 

£ 

k=l 


w {x i - xcg k )} 





1 _ 

2 w, 
k 


(C8) 


(C9) 


Here, the parameters are the local coordinates of the integration points 

with weights w^. The derivatives of the Cartesian coordinates with respect to 
the local coordinate (g,ri) are found from equation (D6) . 
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APPENDIX D 

Representation of the Wing Boundary 

The wing contour c is divided into a series of curvilinear line-elements 
according to Arygris and Scharpf (ref. 17) . A Hermitian interpolation model 
is used to describe the geometry of each curvilinear element. A Lagrangian 
interpolation model is used to describe the vortex sheet intensity over each 
element. Figure D. 1 illustrates a single curvilinear element with the vor- 
tex sheet of intensity y acting over it. Each element is defined by three 
nodes . 



Figure D.l The curvilinear element after Argyris and Scharpf (ref. 17). 


The geometry of a curvilinear element is defined by the known Cartesian 
coordinates of the element end-nodes and the prescribed tangential directions 
at those nodes. The Hermitian polynomials used are 

n i = J ( 2 + 3? - 5 3 ) 

(D1 

H 0 = i (-1 + K + € 2 - £ 3 ) 



(Di) 


H 3 = \ ( 2 - 35 + ? 3 ) 

H 4 = j C *1 * 5 - 5 2 - 5 3 ) 


Here Z, is a. local coordinate defined over the curvilinear element ranging 
in value between -1 and +1. 


A local coordinate point £ is mapped into the (x,y) plane by 

H, 


x(5) 


y(« 



X 1 X I,£ X 3 X 3,? 


y l y l,5 y 3 y 3,? 


H, 


H, 


H, 


(D2) 


The measures of the derivatives of x or y with respect to £ at the nodes of 
the element are given by 


dx 

d? 


k 


x k,£ = q k *x 


dy 

dK 



q . 

i t 

k y 


(D3) 


as suggested by Argyris and Scharpf (ref. 17) . The parameter is defined 

(ref. 17) as 


q 


k 


2 cos <j>k 


(D4) 


and 


cos <f> k = ^ • e 


(D5) 


For an arbitrary point on a curvilinear element, the derivatives of x or y with 
respect to £ are written as 



} 


i 




X l,5 X 3 X 3,S 
y l,5 y 3 y 3,C 






(D6) 
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The derivatives of the Hermitian polynomials with respect to E, are found 
from equations (D2) . 

The intensity of the vortex sheet over a curvilinear element is approxi- 
mated using the Lagrangian interpolation polynomials and the known in- 

tensity of the vortex sheet at the three nodes of the element. At an arbit- 
rary point on the curvilinear element, the intensity of the vortex sheet is 
written as 

3 t e 

yCO = E L.(Oy, = {L}ty} e (D7) 

j=l J J 

where the values of y^ are the nodal values of the vortex sheet intensity. 
The Lagrangian polynomials used in equation (D7) are 

L x = \ ?( 5 * 1 ) 

L 2 = 1 - 5 2 (D8) 

where £ is the local coordinate defined over the element. 
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APPENDIX E 


Representation of the Two-Dimensional Fluid Region 

The two-dimensional fluid region A is represented by subdividing it into a 
finite number of finite-elements. A finite- element mesh is illustrated in fig- 
ure E.l. The flow variables, vorticity w and velocity v=(u,v), are sought 
at the nodes . Two types of finite-elements are used to divide the fluid reg- 
ion A: isoparametric quadrilateral finite-elements and isoparametric triangu- 
lar finite-elements. 



Isoparametric finite-elements better approximate the curved boundary of the 
wing. Considering the local coordinates (£,n) and the corresponding Cartesian 
coordinates (x,y), a point in the plane is mapped into the correspond- 

ing point in the (x,y) plane by 

m 

x (£,n) = E N . (S,n)x (El) 

i-1 1 


and 
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(E2) 


f 


y(5,n) = z N a,n)y 
i=l 

Here, the coordinates (x.,y.) are the known coordinates of the finite- element 

nodes. Each finite-element is defined by m nodes. The shape functions, N^, 

are expressed in terms of the local coordinates which range between un- 

ity and zero within the element. Local derivatives of the shape functions are 
expressed in terms of their Cartesian derivatives by the chain rule. 
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where [J] is the Jacobian. The Cartesian derivatives of the shape functions 
can then be found in terms of their local derivatives as 
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(E5) 


Figure E.2 illustrates the isoparametric quadrilateral finite-element. 
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The shape functions used for quadrilateral finite -element geometry are 

1 2 2 2 2 

N : = j C -l + r - 5n + 5n - 5n + n ) 

n 2 = j ( -l + 5 - 5 n - 5n + +n ) 

N 3 = J ( - 1 + 5 1 2 + 5 2 n + Sn + 5n 2 + n 2 ) 

(E6) 

N 4 = ^ ( -l + 5 2 + £ 2 o - 5n - 5n 2 + n 2 ) 

N 5 = J C l - 5 2 + S 2 n - n ) 

1 2 2 
N, = -=-(1-5 - 5 n + n) 

D Z 

1 2 2 

N ? = j ( l + 5 - 5n - n ) 
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(E6) 


Ng = 2 ( 1 - 5 + 5n 2 - n 2 ) 


Using equations (El) and (E2) , the elements of the Jacobian for the quadrilat- 
eral element may be written as 
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The derivative of the shape functions with respect to the local coordinates is 
computed directly from equation (E6) . 


The isoparametric triangular finite-element is illustrated in figure E.3. 
Local coordinates, used for the triangular finite-element, are the area coordin- 
ates corres P on ding to ^e Cartesian coordinate point (x,y) . The 

applied shape functions are quadratic. The corner nodes are 


N x = £ 1 C 2* - 1 ) 
N 2 " 

N 3 = * 3 C 2* 3 - 1 ) 

and the mid- side nodes 
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(E9) 


The Cartesian derivatives for the triangular finite- element shape functions 
are expressed in terms of their local derivatives. Letting 

5 = n = * 2 1 - ? - n = * 3 (E10) 

then by the chain rule 
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Figure E.3 The isoparametric triangular finite-element. 


Using equations (E10) , the terms of the Jacobian are found from equations 
(Ell) as 
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(E12) 
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The remaining derivatives of the coordinates (x,y) on the right-hand side of 
equation (E12) are found in terms of the area coordinates (£,&,£) using 
equations (El) and (E2) . 12 3 


Over each finite- element, the vorticity and velocity components are approxi- 


mated by 
Then 

their nodal 

values to 

. , u . , and v . 
i i l 

and the nodal 

shape functions 
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m 

= E 
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{N} t {to} e 

(E13) 

and 
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{N} t { v } e 
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The super-script e indicates an element vector or matrix. 
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putational time intervals was assumed to be confined to a thin region around the wing 
surfaces . This report presents thus the foundation and components of the numerical 
method. A following report will present the details of the implementation of the com- 
putational procedure. 
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